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We present an exact dynamical QCD simulation algorithm for the 0(a)-improved Wilson fermion with odd 
number of flavors. Our algorithm is an extension of the non-Hermitian polynomials HMC algorithm proposed 
by Takaishi and de Forcrand previously. In our algorithm, the systematic errors caused by the polynomial 
approximation of the inverse of Dirac operator is removed by a noisy- Metropolis test. For one flavor quark it 
is achieved by taking the square root of the correction matrix explicitly. We test our algorithm for the case of 
Nf = 1 + 1 on a moderately large lattice size (16 3 x 48). The Nf = 2 + 1 case is also investigated. 



1. Introduction 

Lattice QCD simulations with three flavors of 
dynamical quarks are indispensable to under- 
stand the low energy QCD dynamics in the real 
world. Takaishi and de Forcrand [jlj have pro- 
posed a Hybrid Monte Carlo (HMC) algorithm 
which can treat odd number of flavors of dynami- 
cal quarks, in which a non-Hermitian polynomial 
approximation is applied to the inverse of the 
Wilson-Dirac operator. Very recently they have 
removed the systematic error from the polynomial 
approximation, making their algorithm exact 

For realistic simulations with the available 
computational power, the 0(a)-improvement pro- 
gram is widely advocated. We then extend the al- 
gorithm of Refs. to the 0(a)-improved Wil- 
son quark actions. We also develop a noisy- 
Metropolis test to remove the systematic error 
from the polynomial approximation, which is dif- 
ferent from that of Ref. |2j . 

In this article, we describe our algorithm for 
single flavor of dynamical quark. The consis- 

* Presented by K-I. Ishikawa 



tency with the usual algorithm and applicability 
to large-scale simulations are examined by run- 
ning the Nf =1+1 and Nf=2 simulations on a 
16 3 x48 lattice. Employing the algorithm, we 
have started a search for the lattice parameters 
k etc.) suitable for realistic simulations with 
Nf =2+1 flavors of quarks (i.e., degenerate up- 
down quarks and strange quark). We briefly de- 
scribe the novel findings from the search, referring 
to a separate report 0] for details. 

2. Algorithm 

Our HMC algorithm for single flavor of dynam- 
ical quark is based on the following QCD partition 
function: 

Z= [v[U]det{(l + T)]det[D 00 }e- s *M, (1) 

where S g [U] is the plaquette gauge action and 
D oa is the symmetrically preconditioned 0(a)- 
improved Wilson-Dirac operator defined by 

Doo = 1 - (1 + T)-}M oe {l + T)-}M eol (2) 
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where the subscript e (o) means the even (odd) 
components of the Dirac operator. The factor 
T is the local SW-term of the 0(a)-improved 
Wilson-Dirac operator and M is the hopping ma- 
trix. There is another even-odd preconditioning 
in which only the even-even (or odd-odd) SW- 
term is asymmetrically factored out. We com- 
pared the efficiency of the symmetric and asym- 
metric preconditionings within the usual Nf=2 
HMC algorithm, and found that the former has 
better efficiency. Therefore we employ the sym- 
metrically preconditioned form for the QCD par- 
tition function. (For other preconditioning tech- 
niques, see Ref. ||.) 

Using the trick described in Refs. we 
rewrite the determinant of the Dirac operator as 

1 ' d c t[P Npolv [D 00 ]] \dct[T Npoly [D ao ]]\ 2 

where Pn po i v [z] is a non-Hcrmitian polynomial 
of order N po i y (assumed to be even) defined 

by p N foly [z]=Y,i=o'' c i( z - 1 ) i which approximates 
1/z, and W 00 =D 00 P Npoly [D 00 ]. T Npoly [z] is a kind 
of square root of Pn po1v [z] defined by Pn po1 „[z} = 

T* Np Jz}T Npolv [z] withT^JzHE^ 72 di(z-l)\ 
We employ the Chebyshev polynomial approxi- 
mation (i.e., the hopping matrix expansion) for 
Pn po1v [z] and hence the coefficients are Ci={— I) 1 . 
We use the Clenshaw's recurrence formula to cal- 
culate these polynomials. 

Introducing pseudo fermion fields ip a for 
det[T Npoly [D 00 }}\ 2 in Eq. (§), Eq. (@) becomes 

Z = J V[U^ }det[W 00 }e- s ^u My 

S eff [U,ip o ] = S s [U] + S dct [U] + S q [U,ip o ], 

S dct [U] = -log[det[(l + T)]], 

S q [U,*p o ] = \T Npolv [D 00 }^ \ 2 . (4) 

Our PHMC algorithm consists of the following 
two steps in this case: 1) the usual HMC algo- 
rithm for the effective action S e ff, 2) the noisy- 
Metropolis test for the correction term det[Vy oo ]. 
The noisy-Metropolis test is carried out only after 
accepting the HMC Metropolis test. 

The acceptance probability of the noisy- 
Metropolis test is defined by 

P corr {U~+ U'\ =min[f,e- d5 ], 



dS = lAoop'^AooWxol 2 ~ \Xo\ 2 , (5) 

where U is an initial configuration and U' is a trial 
configuration generated by the preceding HMC 
algorithm. Here Xo is a random vector with Gaus- 
sian distribution, and A OQ is defined by A 00 =W 00 - 
This algorithm estimates |detL4 00 ]| 2 instead of 
det [W 00 ] ■ A 00 and A~ Q X are evaluated by the Tay- 
lor expansion with respect to W 00 — 1- In order to 
keep the exactness of our algorithm the residual 
of the Taylor expansion is monitored whenever 
Eq. (||) is calculated. Our programs are written 
in double precision arithmetic, and optimized for 
HITACHI SR8000 at KEK. 

3. N f = 1 + 1 

We test our single-flavor algorithm by ad- 
ditively combining two single-flavored pseudo 
fermions (Nf = 1+1) and comparing results with 
those of the usual Nf=2 HMC algorithm. In this 
case our algorithm runs as follows: 1) the HMC 
algorithm with the two single-flavored pseudo 
fermions, 2) the noisy-Metropolis test for each 
correction terms, which is carried out only when 
the preceding Metropolis test is accepted. To 
check viability of our algorithm toward realis- 
tic simulations, we employ a moderately large- 
scale simulation parameter: 16 3 x48, (3=5.2, 
c sw =2.02, k=0.1340 and k=0.1350. These hop- 
ping parameters correspond to m^/mp ^0.8 and 
~0.7, respectively. 

Figure [I] shows the convergence of T^ poly [D 00 ] 
as Npoiy increases. The measurement is 
made on 20 configurations separated by 10 
trajectories. The residual is defined by 
\ T N poly [Doo}T Npoly [Doo}DooVo-rio\/\Vo\ with a 
Gaussian noise vector r\ - We observe an ex- 
ponential decay as expected and there are no 
accumulation of round-off errors for large N po i y . 

We also investigate the reversibility of the 
molecular dynamics step and observe that the 
violation stays around the limit of double pre- 
cision for both of the quark masses (see Fig. || 
for k= 0.1350). The value of plaquette, averaged 
over ~1000 trajectories, is consistent between 
the two algorithms (e.g., (P)=0.53393(11)(HMC) 
and 0.53392(9)(PHMC with N poly = U0) for 
k=0.1350). With these observation we conclude 
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Figure 1. Convergence T Npoly [D 00 ] as N poly — 
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Figure 2. N po i y dependence of the reversibility 
violations at k = 0.1350. |AiJ|/_ff : total Hamil- 
tonian, AU: gauge link, AP: gauge momentum. 

that our PHMC algorithm for single-flavor dy- 
namical quark works well on moderately large 
lattices and intermediate quark mass for at least 
N f =l+1. 

4. iV> = 2 + 1 

For the Nf =2+1 case, we additively combine 
the standard two-flavored pseudo fermion and 
our single-flavored pseudo fermion as suggested 
in Refs. [^]J^]. The algorithm is given by 1) the 
HMC step with the combined effective Hamilto- 
nian, and 2) the noisy-Metropolis test for the cor- 
rection factor. 

We compare results for the averaged plaque- 
tte between our algorithm and the Hybrid-R al- 
gorithm on a 4 3 x8 lattice at (3=4.8, csw=l-0, 
K M d=0.150, and k s =0.140. Figure || shows 
the molecular dynamics step size dt depen- 
dence of the plaquette value with the Hybrid- 
R and PUMC{N po i y = 10) algorithms. The re- 
sults with the PHMC algorithm (squares) do 
not depend on dt and are consistent with the 
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Figure 3. Plaquette vs dt. 

value at dt — > with the Hybrid-R algorithm 
((P)=0.39669(38)(PHMC), 0.39702(12)(Hybrid- 
R)). 

Encouraged with these results, we haved per- 
formed a series of parameter searches that would 
realize a~ 1 ~l-2 GeV, L~l-2 fm, rn 7r /m p ~0.7- 

0. 8 on a 12 3 x32 lattice. During the search, which 
employed the tadpole- improved one-loop value for 
csw and degenerate quark mass (TV/ =3), we have 
found an unexpected first-order phase transition 
around /3~4. 8-5.0 . This phase transition is also 
observed with the Hybrid-R algorithm on a 8 3 x 16 
lattice with the same lattice parameter. The de- 
tails of the phase structure are given in Ref . || . 
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